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Abstract 

An approximate treatment of exchange in finite-temperature path integral Monte Carlo simula- 
tions for fermions has been proposed. In this method, some of the fine details of density matrix 
due to permutations have been smoothed over or averaged out by using the coarse-grained ap- 
proximation. The practical usefulness of the method is tested for interacting fermions in a three 
dimensional harmonic well. The results show that, the present method not only reduces the sign 
fluctuation of the density matrix, but also avoid the fermion system collapsing into boson system 
at low temperatures. The method is substantiated to be exact when applied to free particles. 

PACS numbers: 02.70.Ss, 31.15.xk, 02.70.-c 
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I. INTRODUCTION 



The path integral Monte Carlo method (PIMC) provides a nonperturbative, basis-set- 
independent, and fully correlated calculation for quantum many-body systems at both zero 
and finite temperature.-^ However for many-fermion systems, PIMC suffers from uncontrol- 
lable errors arising from the notorious sign problem,^ which limits the accuracy or stability 
of the method. The origin of the sign problem comes from the fact that the density matrix 
can be positive or negative by even or odd permutations. At low temperatures, contributions 
from positive and negative parts of the density matrix almost perfectly cancel each other so 
that there is no hope of extracting any useful information. 

A few methods have been proposed to deal with the sign problem. For prob- 
lems in the continuous space, there are fixed-node approximation,- 1 ^ released node 
methods,- exact cancelation methods with Green's function sampling, 8 multilevel block- 
ing algorithm,-*^ a hybrid path integral and basis set method,— various pseudopotential 
approximations,— dlMdkd£.i21J&2M2. the general method for replacing integration over pure 
states by integration over idempotent density matrices,— and method by introducing sev- 
eral images of the system,— global stationary phase approach.— There are also a number 
of methods for lattice models.— So far many efforts have been devoted to tackle this 
problem, however, it remains being the key bottleneck in using PIMC for many-fermion 
systems. 

In this paper, we use a methodology to reduce the rapid oscillation of integrand in the 
evaluation of high dimensional integrals. The idea is that, for the region in which the 
function is rapidly oscillating, the coarse-grained approximations are used to kill fluctuations. 
Applying this technique to PIMC, we found that, the sign fluctuations of density matrix 
can be reduced, thus the Metropolis MC integration algorithm converges efficiently More 
importantly, after carrying out this approximation, the exchange determinant becomes a 
nonlocal form in imaginary time, thus the collapse behavior can be avoided (see below). The 
basic strategy can also be used in the evaluation of other integrals, where the integrands 
exhibit the rapidly oscillating characters. 

This paper is organized as follows. The methodology is described in Sec. II. The numerical 
tests are presented in Sec. III. Discussion and conclusion are given in Sec. IV. 
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II. METHODOLOGY 



To illustrate the coarse-grained approximation used in present paper, we first consider 
the following integral, 

/= / / f(x,y)g(x,y)h(x,y)dxdy. (1) 

J a J A 

We assume that, except for f(x,y) is a rapidly oscillating function for variable y, the rest 
parts of integrands are well behavior (or slow varied) function of x and y. Due to the rapid 
oscillation of f(x,y), it could cause the difficulty on evaluating the integral (Eq. 1) in MC 
simulation. To overcome the difficulty, our strategy is to make coarse-grained approximations 
to f(x,y). To do this, we rewrite the above integral as, 

/= / / F(x)g(x,y)h(x,y)dxdy, (2) 

Ja J A 

with 

F ^ = J A f(x,y')g(x,y')h(x,y')dy' 
!a ' 9(x,y')h(x,y')dy' 

One can see that, F(x) is a kind of coarse-grained functions, and the rapid oscillation of f(x,y) 
is smoothed over. F(x) also can be viewed as an average of f(x,y) weighted by g(x,y)h(x,y). 
If F(x) can be evaluated either exactly or approximately, the rapid fluctuation due to f(x,y) 
could be reduced effectively. For real problems, F(x) usually is hard to be evaluated exactly. 
However, it is possible to determine F(x) under some reasonable approximations, as we have 
done in present paper. 

Considering a three-dimensional system consisting of N spinless, indistinguishable quan- 
tum fermions, the standard PIMC is based on the following expansion of partition function: 

„ N M 



M mZo I II fl dfj v) detAexp(-pH) (4) 
~*°° J i=lv= i 



with 



N , M 



i=l v=\ 

where ujm = ^M't C = )'^ NM , r^ M+1 ^ = f^ 1 -*, V is the potential energy, and the 

square length(Lj) is defined as, Lf = Yl^i^/f^ ~ The subscript i refers to the 

particle number while the superscript v refers to different slit of imaginary time. (3, m 
and M are reciprocal temperature (l/ksT), mass of particles and total number of beads 



3 



respectively, {rf } refers to (r^\r£\ ,f^_ v ry).To make the expression compact, we 

have introduced a N x N matrix A whose element reads 



Ay = exp{-\(5mu>\M > _ r f y _ {r f) _ r f y )y (5) 



detA is the determinant of the matrix A, which accounts the contribution of permutations to 
the partition function. It is detA, which can be positive or negative, that causes the so-called 
sign problem. Previous studies based on pseudopotential methods have shown that, direct 
use of Eq.4-like formula usually results in a fermion system collapsing into a bosonic state at 
low temperature."^ The physical reason comes from the fact that the matrix A approaching to 
unit at low temperature. To prevent this undesirable behavior, people usually recast matrix 
A in a nonlocal form as suggested by Hall, 14 or directly use a nonlocal pseudopotential as 
suggested by Miura and Okazaki.— Although these schemes do give a good solution, the 
computational cost also increases. 

From Eq. 5, one can see that, if r\ M ^ is close to ij\ Ay will be a rapidly oscillating 
function of rp ■ Although, for N>2, it is difficult to prove the direct relation between the 
rapid oscillation of Ay and the sign problem, it is quite clear that the rapid oscillation of /In- 
directly results in the sign fluctuation for N=2. To smooth the rapid oscillation, the coarse- 
grained approximation is made for Ay by integration over for all possible configurations 
with fixed r\ M \ rf 1 ^ and Lj. Now Eq. 4 is replaced by 

„ „ N M 

Z = W\ ifeo / II II M U) detEexp(-(3H) (6) 

1=1 v=\ 

where the new N x N matrix H is the coarse-grained approximation of matrix A. According 
to the idea presented in Eq. 1, 2 and 3, the elements of H are 

f'Ut=i 1 drf ) A a ,exp(-pH) 



I'UtL^drVexpi-PH) 

Since A ai is only relevant to and fj , Eq. 7 can be rewritten as, 

„ _ f U.t=i drf ] A ai exp(-\f3mul 4 Ll - (3U) 
^ f dr^expi-lPmulLl - (3U) 



(7) 



(8) 



with 



M 
u=l 



J' in Eq. 7 and 8 refers the integral under the constraint of fixed r£^\ f^ M ^ and L a . 

Since the kinetic energy relevant part^cc^L 2 ) is a constant for fixed L a . The Eq. 8 can 
be further simplified as, 

- _ /' nf^i 1 dfi u) A a ,exp(-pU) 



-•cry 



I'U^dr^expi-pU) 
To calculate the Eq. 9 under the constraint of fixed L a , we can rewrite it as, 



(9) 



(10) 



where A(L a , r£\ ri^) = /' FI^ 1 dr£\ which is the number of configurations for fixed L a , 
and . U reads, 

. M-l 

U= / ]ldffexp(-pU)/A(L a ,i$\rM) (11) 

The above integral is also under the constraint of fixed L a . There is almost no hope to 
evaluate Eq. 10 exactly. In this paper, we calculate it with approximations, 

2 a7 « — tA(m)7 / d^A ai A(L a , t£>, r^), (12) 
T(L a ,ri M ^) is the total number of configurations for fixed L a and ri M \ i.e, 



r(L a ,fi M) ) = jdr^A(L a ,rt\r^). 



By replacing Eq. 10 with Eq. 12, we have assumed that U is weakly dependent of 
for fixed ri M \ r) M ^ and L a . This approximation works well if M is not too small. The 
reason lies on the fact that, only the configurations, in which is close to ri M \ make the 
A(Lj, rf 1 ^ 1 ) be significant (see below Eq. 13 and 14). Our numerical test (see below) 
also demonstrates this point. 

A(L a , r£\ ri M ^) can be written as an integral over three Cartesian directions, 

A(WlVl M) ) = J A(L a ^),xW)A(L a ^ 

(13) 

where the integral is evaluated under the constraint: L 2 a = L 2 ax +L 2 ay +L 2 az . A(L ax , Xa\x^) 
is the number of configurations for fixed L ax (L 2 ax = J2^f =1 (xa +1 ^ — x^) 2 ) ,ia' and x^\ 
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A(L ayj ya\y^) and A(L az , z£\ z a ) are the counterparts of A(L ax , Xa\ x^) along y 
and z direction respectively. 

To calculate A(L ax ,Xa ,x t), we define a (M-f )-dimensional vector R, of which Carte- 
sian components are (xa — %a \ Xa 4 ^ — ii M '). First, for a given L ax , it requires 
\R\ = \ I? — (xa — Xa ) 2 i all the configurations satisfying this condition lie on a surface of 



(M-l)-dimensional super-sphere with radius equal to y L^ x — (x$ — x^) 2 ; Second, since 
the Cartesian components of R are not independent, i.e, the project of R on the (M-l)- 

/ (1)_ (M) \ 

dimensional unit vector is , this condition defines a (M-l)-dimensional super-plane; 

Thus, all the (M-l)-dimensional points, which attribute to A(L ax , Xa\x^), lie on a (M-2)- 
dimensional super-spherical surface intersected by the (M-l)-dimensional super-sphere and 
the (M-l)-dimensional super-plane. According to analytic geometry in high dimensional 
space, A(L aX} Xa \ x^) is the proportional area of the (M-2)-dimensional super-spherical 

/ (1)_ (M)n 2 1 

surface with radius equal to (R 2 — ) 5 - We end up with: 

/ (1) _ (M)x 2 

A(L ax ,xW,xM)dRdxU oc C A (R 2 - J ^dRdx® 

with = (M — 2)ti M 2 1 /((M — 2)/2)!.— By changing integration variable from (R, x l a ) to 
(L ax , x a ), we have, 



A(L ax ,x^,x^)dL ax dx^ oc C A L ax (L 2 ax - (x« - x^ )f dL ax dx^\ (14) 

Similarly, we can obtain A(L ay , yi , yk ) and A(L az , Za , z Q M ), which have the same 
formula as A(L ax , Xa , Xa )• Substituting A(L ax , io',i[, M ') in Eq. 14 for A(L ax , x& ,Xa) 
in Eq. 13, as well as replacing the counterparts of A(L ay ,ya\y^) and A(L az , Za \ z a ) 
in Eq. 13, A (L a , r£\ r£ ) can be obtained. As a result, E ai ought to be calculated 
numerically. 

One can see that, A(L ax , Xa ) quickly decays as a function of \x& — x a ; |. The be- 
haviors of A (L^, , 2/a ) and A(L a2 , z^, Zcf^) are the same as that of A(L ax , Xa \ Xa^). 
Accordingly, A(L a , , fa^) also quickly decays as a function of \fa — rl M ^|- By employ- 
ing the change of variables as ri 1 ^ = rl M ' ) + 5 a , and making coordinate transformation in 
spherical coordinates, we can see that S a7 is a function of L a and r&P (r^P = \fc^ — r« |). 

After carrying out above coarse-grained approximation, we find that, the off-diagonal 
element H a7 is a function of both rffl and L a , which has an explicit nonlocal form. In 
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contrast, the off-diagonal element A ai is a local function in imaginary time. More impor- 
tantly, the off-diagonal element A ai could be much larger or much smaller than 1.0, while 
the off-diagonal element S Q7 is less than 1.0 when L a is not too long (see Fig. 1). Thus, 
after replacing A ay with S a7 , at least for two-particle system, a lot of sign fluctuations are 
well canceled. However, since the length of path could be much longer at low temperature, 
the off-diagonal element S a7 can be larger than 1.0 (see Fig. 1), which will result in the 
presence of negative detE. To further reduce the negative sign, we have made the second 
stage of coarse-grained approximation for longer paths. Using the similar idea as above, we 
can replace Eby & N x N matrix A, 



A Q7 = { B r^° rL " ~ K (15) 



L* is a function of which is determined by, 

The above integral is made over all the configurations with L a > L* a . For an arbitrary 
interact potential, the above integral is almost no hope to be evaluated exactly. However it 
can be calculated with approximations, 

^ dL a E ai exp(-^(3mw 2 M L 2 a )T(L a , fj M) ) 

In the evaluation of the above coarse-grained approximation, we have assumed that total 
potential energy is constant when L a is longer than certain value, i.e, L* a . In current work, 
we have taken rj(r^) = 1 through the whole paper. 

Within the current approximations, A has a few advantages over the original matrix A. 
First, our calculations have shown that the off-diagonal element of A is not larger than 
1.0 anywhere (see Fig. 1), in constrast the off-diagonal element of A could be much larger 
than 1.0. Thus at least for two-particle system, A is always non-negative, the sign problem 
completely vanishes. Second, different from A, A is nonlocal, which depends on the whole 
path. The nonlocal behavior of A can effectively avoid the collapse of fermion system into 
boson system at low temperature. At low temperature, the length of path becomes longer 



and longer, so the off-diagonal element of A has more chance being 1.0. This situation 
makes A have little chance being unit matrix. Our calculation also demonstrates this point. 
It should be pointed out that the current formula is exact for free particles (see APPENDIX). 
Now we end up with the final formula for real calculations, 

„ N M N M 1 M 

z = ° / n n w e e - ^ ) b + e 

i=l u=l i=l v=\ u=l 

In real calculation, the element of A is first numerically integrated. At the same time, 
the derivative of A respective to temperature is also numerically calculated to account the 
contribution to thermal energy. Eq. 18 can not be directly used in standard MC, since for 
the fermionic systems, detA is not always positive. However Eq. 18 can be integrated using 
modified MC technique, which is widely used previously.— 1 ^ To achieve this result, we first 
defined pseudo-Hamiltonian, H p , which is, 

N M 1 M 

E v = E E 2™m^ +1) - *f V + E ^m^}) + ln\detA\. 

i=l u=l v=\ 

The thermodynamic average of a physical quantity Q is 

(Q) = I Uli UtU Dfj u) Q({r ( i , ' ) })sgn(detA)exp(-(3H p }) m 
I nil UtU D^sgn(detA)exp(-(3H p })}) ' 

where sgn(detA) stands for the sign of detA at a configuration. 

It needs to be pointed out that, we have used the similar technique as most pseudopo- 
tential methods,— 1 i2 1 Mii^iiSiiLi§ J iSi2S but we do not recast matrix A or extend A into each 
imaginary time. 



III. THE NUMERICAL TESTS 



To illustrate the usefulness of the current method, we have considered N interacting 
spinless fermions confined in a three-dimensional harmonic well, which Hamiltonian reads, 

N si 2 N 

^B^ + ^D+Enr,), (20) 

j=l i<j 

where m, fj, pj and V(ry) are mass, positions, momenta of the particles, and the inter- 
particle interaction potential, respectively. For computational simplicity, the units by which 
m = h = ks = 1 are used through the rest of this paper. In current calculations, we 
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consider three cases, i.e, Case 1: V(ry) = 0, N=6 and uo 2 = 1, no interaction between 
particles, reflecting a standard harmonic system; Case 2: N=6, Vfcj) = —^j-rfj, where 
the interaction is also harmonic one with Q 2 = and u 2 = 4; Case 3: V(ry) = 
uo 2 = 0.320224986, q = 1 and N=2, interaction between particles is the Coulomb potential, 
the parameters correspond to hydrogen-like ion (H + ) of Kestner-Sinanoglu model. 31 The 
exact results of all three cases can be found elsewhere,— 1 ^ which is easy to check the validity 
of the current method. These models are widely used as a benchmark for checking the 
usefulness of various methods for sign problem, see for examples Ref.— i22i2ii22, 

Our Metropolis MC scheme is preformed based on Eq. 19. At each step, H p are calculated 
to determine the rejection and acceptance. detA is calculated by a certain algorithm with 
the computational cost scaled by N 3 . 28 This kind of numerical technique enables us to 
perform the fermionic simulations with reasonable computational time. There are two basic 
types of moves in current simulations: (1) Displacement move, where all the coordinates 
for a single particle are displaced uniformly; (2) Standard bisection moves.-^ 3 - The MC 
procedure used in this work is wildly used by others. One MC step is defined as one 
application of each procedure. Ten million MC steps of calculation were carried out for 
each temperature. For a few cases, 100 million MC steps are made to check the ergodic 
problem. The results agree with the short runs within the error bars. To further check 
the ergodic problem, the simulations are carried out by a few random generated starting 
configurations. All simulations converge to the same results. The energy is calculated based 
on the thermodynamic estimator.- The energies are well converged at M//3=20, 22 and 5 
for case 1, 2 and 3 respectively. 

The calculated thermal energy is in good agreement with the exact one for all three cases 
studied. In case 3, the exact energy 2.647 of Ref.— has been almost accurately reproduced, 
which is 2.652±0.003 in current simulations. Fig. 2 shows the thermal energy per particle as 
a function of temperature for Case 1 and 2, the corresponding exact results are also shown 
in Fig. 2 with lines. As can be seen from the figure, the overall temperature dependence 
is well reproduced by current calculations. The calculated thermal energies agree very well 
with the exact value at low temperature. The slight deviation at high temperature is due to 
the fact that the first stage of approximation will result in error when the number of beads 
is too small, which is the case for high temperature. 

We have calculated the pair correlation function (PCF) between beads, which is defined 
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2 



M N-l N 



MN(N - 1) 



E££*( 




v i j>i 



It is known that,— comparing with boson and Boltzmann systems, the fermionic PCF has 
a hole around the origin, which reflects the Pauli exclusion principle. In Fig. 3, we present 
PCF for case 1 and 2 at temperature of 0.2. From this figure, we can see that, the pair 
correlation function clearly represents the effect due to the Pauli exclusion principle. The 
similar behaviors are observed for other temperature and systems. 

The average sign reflects the signal-to-noise ratio, which directly affects calculation preci- 
sion and computation time needed. The average sign is defined as Sign = (N + — AL)/ (N + + 
N-), where iV + and AL are the total positive and negative configuration respectively. The 
lower panel of Fig. 4 shows the average sign of current simulations via temperature. It can 
be seen that, Sign decreases with the decrease of the temperature. However, for the studied 
systems, even at lowest temperature (T=0.1), the average sign is quite high (around 0.1). 
We also calculated the average sign via the number of particles at T=0.5 for both case 1 
and 2, which is shown in the upper panel of Fig. 4. Similarly, Sign also decreases with 
the increase of the number of particles. Although we have not completely solved the sign 
problem, our approach does much improve the sign decay rate with both temperature and 
number of particles. Direct using of Eq. 4, Sign is about 0.01 at temperature of 0.8 for 
case 1. And for temperature lower than 0.8, the large sign fluctuation makes MC simulation 
difficult to obtain any useful information. According to the data shown in Fig. 4, the max- 
imum number of particles, which can be handled in current method, should be in order of 
ten. Considering both spin-up and -down, the maximum number of particles can be around 
twenty, which could be particularly useful for atom and molecular systems. 

IV. DISCUSSION AND CONCLUSION 

In this paper, we have introduced an approach to reduce the fermion sign fluctuation in 
finite temperature PIMC simulations. By this method, configurations, which probably cause 
the sign fluctuation, are pre-calculated within two stages of coarse-grained approximations, 
while the rest are treated exactly. After two stages of coarse-grained approximations, at 
least for two-particle system, the sign problem is solved completely. Since the exchange 
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matrix A is replaced by a non-local one (A), the collapse of fermion system into a boson 
one at low temperature has been effectively avoided. The pilot calculation was performed 
on three model systems: six independent particles in a three-dimensional harmonic well, six 
interacting particles in a three-dimensional harmonic well, and hydrogen-like ion (H + ) of 
Kestner-Sinanoglu model. The calculation shows that the current approach not only dra- 
matically drops the sign fluctuation, but also gives an excellent description to real systems. 
Our method could be particularly useful for atom and molecular systems. Although our ap- 
proach suffers from the sign problem for large number of particles, we believe that it provide 
an alternative thought on the sign problem. We also believe that a similar approach can 
also be helpful in other path integral methods. The current formula can be easily extended 
to systems consisting of both spin- up and -down fermions. (see for example,— 1 ^) 

Our approximation breaks down for systems including particles more than twenty (in- 
cluding both spin up and down particles). It would be possible to generalize our method for 
problems of larger numbers of fermions. Although we have used ^{r^ 1 ) = 1 through out 
this paper, other values are also possible. For example, if T)(r^) = exp{—a^muj\ I (r\^) 2 ) 
is chosen, the current method can be more flexible. For a = 0, it is the case used in current 
work. For a — 1, detA becomes the exact density matrix of free particles(see APPENDIX), 
thus the sign problem can be avoided completely. In fact, with a increasing from to 1, the 
approximation becomes more and more crude, but the negative parts become less and less. 
To further improve the current method, a better form or value for rj(r^) could be found. 
It is actually the issue on which we are working now. 
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APPENDIX 

In this appendix, we will prove that the current formula is exact for free particles. Since 
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the second stage of approximation is just a straightforward integration for free particles, we 
only prove the formula of first stage is correct for free particles. All Cartesian coordinates 
are equivalent for free particles, for simplicity we only prove it in one Cartesian direction, 
say, x. 

The partition function for free particles in one dimension has the form, 

z=^/n^<%({4 M, }.{4 M) }) (2D 

,W (M) 
J i i x i 

P»({x?h {*?}) = Cm [ I] IT Dx^ D) exp(-P £ \rmJ M L\ x ) (22) 

J i=l v =\ i=l 

where C\d = ( 2 ^^ 2 )^, and 2^ D ' ) is the one-dimensional counterpart of S^-, which is 



i=l 

where p(x\ M> ,x\ M> ) is the density matrix, of which element with current formula reads, 

N M-l N 



"13 TIT ( M )~\ 



(23) 



T 1D (L fa ,x^) = j dx?K{L lx ,xf\x^) 

Eq. 22 is only relevant to {L ix } and {x^}, the integration over x^ (z/=l,...M-l) can be 
replaced by (M-l)-dimensional spherical polar coordinates, i.e, integration over L ix multi- 
plying T 1D (L ix ,xf), Eq. 22 becomes, 

Puds?}, K M » « C 1D J dL xl dxfK{U x , xfKx^) (24) 



xexp(--(3mcol(Ll + (^ - ^ " & ' - 
Remembering for fixed x^p and x[ M \ the minimum value of Lf x equals to (x^ — 
x i M ^) 2 WZT- We first do the integration over variable {L ix } by change of variables as 
L% = Ll - (*S 1} - xf) 2 ^, Plj ({xf }, {xf}) becomes, 



A,({*f }, R M }) oc g lg C A ^ ^ y ^ (25) 



■ M-4 
■ 2 

xexp(-i/3m^((x« - xf } ) 2 + - z™)*)) 



gl ^ ( ^ } ^ V M(3muj 2 M ^P(~2M^M ~A ) 
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where Cm is an irrelevant constant. The above express is the exact formula for free particles. 
It needs to noted that, since we only can get the relative value for A, we could not obtain 
the absolute value of Cm- However the absolute value of Cm is irrelevant to our calculation, 
which is a function of M only. Our numerical test also shows that the calculated element of 
density matrix based on the current formula is in excellent agreement with the exact data. 
In fact, the current formula must be exact for free particles, since all the approximations 
become exact without potential part. 
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1.2 




FIG. 1: (color online) The off-diagonal element of A and S (Ajj and Sj,-) as a function of particle 
separation (rij) for a few selected lengths of path (L) at T=0.5. When the length of path is short, 
Ajj and Hjj are the same. As length being longer, Ejj can be larger than one for short particle 
separation. In contrast, Ajj is not larger than one for any case. 
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FIG. 2: (color online) The thermal energy of case 1 (Circle) and 2 (Squares) as a function of 
temperature calculated by PIMC. The solid and dash lines are the exact energies of case 1 and 2 
respectively. The agreement is quite well. 
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FIG. 3: The pair correlation functions between beads at T=0.2 for both case 1 and 2. The hole 
around original is the reflection of Pauli exclusion principle. 
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FIG. 4: The average sign of case 1 (circles) and 2 (squares) via temperature (lower panel) and 
number of particles (upper panel) as calculated by our PIMC simulation. Sign decreases with the 
decrease of temperature, and the increase of number of particles. 
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